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ABSTRACT 

Wc present 2-D hydro-dynamic simulation of rotating galactic winds driven by radia- 
tion. We study the structure and dynamics of the cool and/or warm component (T ~ 
10 4 K) which is mixed with dust. We have taken into account the total gravity of a 
galactic system that consists of a disc, a bulge and a dark matter halo. We find that 
the combined effect of gravity and radiation pressure from a realistic disc drives the 
gas away to a distance of ~ 5 kpc in ~ 37 Myr for typical galactic parameters. The 
outflow speed increases rapidly with the disc Eddington parameter ro(= k//(2c(7S)) 
for To > 1.5. Wc find that the rotation speed of the outflowing gas is < 100 km s _1 . 
The wind is confined in a cone which mostly consist of low angular momentum gas 
lifted from the central region. 
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1 INTRODUCTION 

Many galaxies are observed to have moving extraplanar gas, 
generally termed as galactic superwinds (see Veilleux et al. 
2005 for a recent review). Initial observations showed the Ha 
emitting gas above the plane of M82 (e.g. Lynds & Sandage 
1963). The advent of X-ray astronomy established yet an- 
other phase of galactic outflows, namely the hot plasma, 
emitting X-rays in the temperature range 0.3-2 keV (Strick- 
land et al. 2004). Also recent observations have revealed the 
existence of molecular gas in these outflows (Veilleux et al. 
2009, waiter et al. 2002). Earlier observations were limited 
to local dwarf starburst galaxies that showed these winds. 
However, in recent years, the observations of outflows in Ul- 
tra Luminous Infra-red Galaxies (ULIGs) have extended the 
range of galaxies in which outflows are found (Martin 2005, 
Rupke et al. 2005, Rupke et al. 2002). 

On the theoretical side, there have been speculations 
on winds from starburst galaxies (Burke 1968, Mathews & 
Baker 1971, Johnson & Axford 1971). In these models the 
large scale winds are a consequence of energy injection by 
multiple supernovae (Larson 1974, Chevalier & Clegg 1985, 
Dekel & Silk 1986, Heckman 2002). In the context of the 
multiphase structure of the outflows, the results of these 
theoretical models are more relevant for the X-ray emitting 
hot wind. On the other hand, observations of the cold out- 
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flows are better explaind by the radiation driving (Murray 
et al. 2005, Martin 2005). 

If only Thompson scattering is considered, then radia- 
tion from galaxies does not seem to be a reasonable wind 
driving candidate because opacities would be small; how- 
ever one should consider that these winds are heavily en- 
riched. Murray et al. 2005 proposed a wind driving mecha- 
nism based on the scattering of dust-grains by the photons 
from the galaxy (see also Chiao & Wickramasinghe 1972; 
Davies et al. 1998). This mechanism can be quite effective 
since the opacities in dust-photon scattering can be of the 
order of hundred cm 2 g _1 and gas in turn, being coupled with 
the dust, is driven out of the galaxy if the galaxy posseses a 
certain critical luminosity. Bianchi & Ferrara (2005) argued 
that dust grains ejected from galaxies by radiation pressure 
can enrich the intergalactic medium. Nath & Silk (2009) 
then described a model of outflows with radiation and ther- 
mal pressure, in the context of outflows from Lyman break 
galaxies observed by Shapely et al. (2005). Murray et al. 
(2010) have also described a similar model in which radia- 
tion pressure is important for the first few million years of 
the starburst phase, after which SN heated hot gas pushes 
the outflowing material. Sharma & Nath (2011) have also 
shown that radiation pressure is important for outflows from 
high mass galaxies with a large SFR (with v c 200 km s~ , 
SFR ^ 100 M yr" 1 ), particularly in ULIGs. 

In this paper, we study the effect of radiation pressure 
in driving cold and/or warm gas outflows from disc galaxies 
with numerical simulations. Recently, Sharma et al. (2011) 
calculated the terminal speed of such a flow along the pole of 
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a disc galaxy, taking into account the gravity of disc, stellar 
bulge and dark matter halo. They determined the minimum 
luminosity (or, equivalently, the maximum mass-to-light ra- 
tio of the disc) to drive a wind, and also showed that the 
terminal speed lies in the range of 2-4 V c (where V c is the ro- 
tation speed of the disc galaxy) , consistent with observations 
(Rupke et al. 2005, Martin 2005), and the ansatz used by 
numerical simulations in order to explain the metal enrich- 
ment of the IGM (Oppenheimer et al. 2006). We investigate 
further the physical processes for a radiation driven wind. 
Rotation is yet another aspect of the winds that we address 
in our simulation. As the wind material is lifted from a ro- 
tating disc, it should be rotating inherently which is seen in 
observations as well (Greve 2004, Westmoquette et al. 2009, 
Sofue et al. 1992, Seaquist & Clark 2001, Walter et al. 2002). 

Previous simulations of galactic outflows have consid- 
ered the driving force of a hot ISM energized by the effects 
of supernovae (Kohji & Ikeuchi 1988; Tomisaka & Bregman 
1993; Mac Low & Ferrara 1999; Suchkov et al. 1994, 1996 
; Strickland & Stevens 2000; Fragile et al. 2004; Cooper et 
al. 2008, Fujita et al. 2009). However the detailed physics 
of a radiatively driven galactic outflow is yet to be studied 
with a simulation. In this work, we study the dynamics of an 
irradiated gas above an axisymmetric disc galaxy by using 
hydrodynamical simulation. Recently Hopkins et al. (2011) 
have explored the relative roles of radiation and supernovae 
heating in galactic outflows, and studied the feedback on 
the star formation history of the galaxy. Our goal here is 
different in the sense that we focus on the structure and 
dynamics, particularly the effect of rotation, of the wind. In 
order to disentangle the effects of various processes involved, 
we intentionally keep the physical model simple. For exam- 
ple, we begin with a constant density and surface brightness 
disk, then study the effect of a radial density and radiation 
profile, and finally introduce rotation of the disk, in order 
to understand the effect of each detail separately, instead 
of performing one single simulation with many details put 
together. 



2 GRAVITATIONAL AND RADIATION 
FIELDS 

The main driving force is radiation force and the containing 
force is due to gravity. We take the system to be composed 
of three components disc, bulge & dark matter halo. We 
describe the forces due to these three constituents below. 
We take a thin galactic disc and a spherical bulge. All these 
forces are given in cylindrical coordinates because we solve 
the fluid equations in cylindrical geometry. 



2.1 Gravitational field from the disc 

Consider a thin axisymmetric disc in r<f) plane with surface 
mass density E(r). As derived in the Appendix, the vertical 
and radial components of gravity due to the disc material at 
a point Q above the disc with coordinates (r, 0, z), are given 
by 
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Figure 1. Magnitude of gravitational force of the (a) uniform 
disc (UD) (b) exponential disc (ED) in colours with direction in 
arrows. Values are in the units of GSq(= 4.5 X 10 -9 ) dyne. 
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(1) 



The azimuthal coordinate of Q is taken to be zero, because 
of axisymmetry. The integration limit for cj>i — to 2-7T. 

We consider two types of disc in our simulations, one 
with uniform surface mass density and radius rd (UD), and 
another with an exponential distribution of surface mass 
density (ED) with a scale radius r s . The surface mass density 
of uniform surface density disc (i. e.,UD) is 



S = Eo = constant 

and in the case of a 
distribution (ED) 

E = E exp(-r//r s ), 



(2) 

disc with exponentially falling density 



scale length . 



(3) 



In case of UD (eqn [2} , the integration limit would be rl — 
to rd, while for ED (eqn [3}, the limits of the integration 
run from rl — to oo. Numerically this means, we integrate 
up to a large number, increasing which will not change the 
gravitational field by any significant amount. We have cho- 
sen the Es in such a way that the total disc mass remains 
same for the UD or ED. Therefore, 



(4) 



+ z 2 + rl 2 



2rrlcos<j>l] 3 / 2 



In Figure 1, we plot the contours of gravitational field 
strength and its direction vectors due to a UD (left panel), 
and that for the ED (right panel). Interestingly, discs with 
same mass but different surface density distributions, pro- 
duces different gravitational fields. For the UD the gravita- 
tional field is not spherical and the gravitational acceleration 
is maximum at the edge of the disc. On the other hand, the 
field due to ED is closer to spherical configuration with the 
maximum being closer to the centre of the disc and falling 
off outwards. 



2.2 Bulge and the dark matter halo 

We consider a bulge with a spherical mass distribution and 
constant density, with mass Mb and radius r^- The radiation 
force due to the bulge is negligible as it mostly hosts the old 
stars. The gravitational force of the bulge is given by 
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otherwise 



GM h z -c r, , 
—J 1 - , lfR< r b 



- G t I .t" , otherwise 



(5) 



(6) 



where R = y/r 2 + z 2 . 

We consider a NFW halo with a scaling with disc mass 
as given by Mo, Mao and White (1998; hereafter referred to 
as MMW98) where the total halo mass is ~ 20 times the 
total disc mass. The mass of an NFW halo has the following 
functional dependence on R 



M(R) = 4m-poritS R s 



where x — „ R 

-K200 



ln(l + cx) 



1 + cx 



-^200 



So = 



(7) 



Here 



3 !n(l + c)-c/(l + c) ■ 

pcrit is the critical density of the universe at present epoch, 
R s is scale radius of NFW halo and R200 is the limiting 
radius of virialized halo within which the average density is 
200 pcrit - This mass distribution corresponds to the following 
potential, 



&NFW — —4.TY pcritSoRg 



In (1 + R/R s )/R 



(8) 



The gravitational force due to the dark matter halo is there- 
fore given by, 



fhalo,r 
fhalo,z 



dr 

d&NFW 



r GM(R) 

(r 2 + z 2 y/ 2 ' 

z GM(R) 



dz (r 2 + z 2 ) 3 / 2 ' 

The net gravitational acceleration is therefore given by 
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The gravitational field for both bulge and halo is spher- 
ical in nature, although, that due to the bulge maximises 
at r-j,. However, the net gravitational field will depend on 
the relative strength of the three components. In Figure [2] 
(left panel), we plot the contours of total gravitational field 
strength due to the bulge, the halo and an UD. The non- 
spherical nature of the gravitational field is evident. A more 
interesting feature appears due to the bulge gravity. The 
net gravitational intensity maximizes in a spherical shell of 
radius r^(= 0.2L re f, see section §3.1). Therefore, there is 
a possibility of piling up of outflowing matter at around a 
height z ~ r b near the axis. In the right panel of Figure j2]), 
we present the contours of net gravitational field due to an 
embedded exponential disc within a halo and a bulge. 



2.3 Radiation from disc and the Eddington factor 

We treat the force due to radiation pressure as it interacts 
with charged dust particles that are assumed to be strongly 
coupled to gas by Coulomb interactions and which drags the 
gas with it. The strength of the interaction is parameterized 
by the dust opacity k which has the units cm 2 gm _1 . 

Gravitational pull on the field point Q(R, Z) due to 
the disc point P(rf, <j)f, 0) is along the direction (sec 
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Figure 2. Total gravitational force of the (a) uniform disc (b) 
exponential disc in colors with direction in arrows. The values are 
in the same units as in Figure [T] 



appendix). The difference in computing the radiation force 
arises due to the fact that one needs to account for the pro- 
jection of the intensity at Q (for radiation force from more 
complicated disc, see Chattopadhyay 2005). For a disc with 
surface brightness I(r), we can find the radiation force by 
replacing GE(r/) in eqn[T]by I(rt)n/c, and take into account 
the projection factor z/^/r 2 + z 2 + rl 2 — 2rrf cos <j)f. Similar 
to the disc gravity, the net radiation force ~^ r ad at any point 
will have the radial component (F ra d, r ) and the axial com- 
ponent (F ra d, z ) and are given by, 

d(/>fdrtl(rt)(r — rtcostf)/) rl 



F, 



r(r,z) 



ffl 



r 2 + z 2 + rl 2 — 2rr/cos(f>/] 2 



F ra d,z(r,z) 



KZ 
C 



fr,r(r,z) 



d(j>ldrll{rl)rl 



[r 2 + z 2 + rl 2 — 2rr/cos(f>/] 



(11) 



(12) 



^rfr,z(r,z) 



Since we have two models for disc gravity, we also consider 
two forms of disc surface brightness. 



I = Io = constant, for UD 
and 

/ = f exp(-r//r s ) ,for ED 



(13) 



(14) 



If the two disc types are to be compared for identical lumi- 
nosity, then one finds 

-\ ■ (15) 



The disc Eddington factor is defined as the ratio of 
the radiation force and the gravitational force (MQT05). 
In spherical geometry this factor is generally constant at 
each point because both gravity and radiation has an in- 
verse square dependence on distance. Although in the case 
of a disc, the two forces have different behaviour, we can 
still define an Eddington parameter as F = p rad . In this 
case this parameter depends on the coordinates r, <f>, z of the 
position under consideration. We can however define a pa- 
rameter whose value is the Eddington factor at the centre 
of the disc, i.e., 

kI 



To = 



2cG£ 



(16) 
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Figure 3. Magnitude of force due to radiation from the (a) 
uniform disc, (b) exponential disc for To = 0.5, with arrows for 
direction. 



If To = 1, then the radiation and gravity of the disc will can- 
cel each other at the centre of the disc. We will parameterize 
our results in terms of To. Therefore, the components of the 
net external force due to gravity and radiation is given by 



T^t- — Fg ra v ,r F ra d,r — GSq \fg,r fr,r) 

Ti-z — F gr av,z — F ra d >z = GEq (fg, z — 2ro/r,z) 



(17) 



In Figure [3] we plot the contours of radiative acceler- 
ation from an UD, and the same from an ED. There is a 
significant difference between the radiation field above an 
ED and that above an UD. While the radiation field from 
an UD is largely vertical for small radii, but starts to diverge 
at the disc edge, at r ~ rd- One can therefore expect that 
for high enough I, the wind trajectory will diverge. In case 
of ED, the radiation field above the inner portion of the disc 
is strong and decreases rapidly towards the outer disc. 



3 NUMERICAL METHOD 

The hydrodynamic equations have been solved in this paper 
by using the TVD (i. e., Total Variation Diminishing) code, 
which has been quite exhaustively used in cosmological and 
accretion disc simulations (see, Ryu et al. 1993, Kang et al. 
1994, Ryu et al. 1995, Molteni et al. 1996) and is based on 
a scheme originally developed by Harten (1983). We have 
solved the equations in cylindrical geometry in view of the 
axial symmetry of the problem. This code is based on an 
explicit, second order accurate scheme, and is obtained by 
first modifying the flux function and then applying a non- 
oscillatory first order accurate scheme to obtain a resulting 
second order accuracy (see, Harten 1983 and Ryu et al. 1993 
for details). 

The equations of motion which are being solved numer- 
ically in the non-dimensional form is given by 



9q 



1 <9(rFi 
r 



dF 2 + dG 

dr dr dz 



where, the state vector is 
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(20) 



and the source function is given by 


r 

-pTZ z 
_ -p[v r TZ r + v z TZ z )_ 



(21) 



3.1 Initial and boundary conditions 

We do not include the disc in our simulations and only con- 
sider the effect of disc radiation and total gravity on the 
gas being injected from the disc. We choose the disc mass 
to be Md = 10 11 Mq and assume it to be the unit of mass 
(i. e.,M re f). The unit of length (i. e.,L re f) and velocity (i. 
e.,v re f) are Td = 10 kpc and v c = 200 km s _1 , respectively. 
Therefore, the unit of time is t re f — 48.8Myr. We intro- 
duce a normalization parameter f; such that GMd/v^ = £7^, 
which turns out to be £ = 1.08. Hence the unit of density 
is p ref = 6.77xl0" 24 g cm -3 (~ 4m p cm -3 ). All the flow 
variables have been made non-dimensional by the choice of 
unit system mentioned above. 

It is important to choose an appropriate initial condi- 
tion to study the relevant physical phenomenon. We note 
that previous simulations of galactic outflows have consid- 
ered a variety of gravitational potential and initial ISM con- 
figurations. For example, Cooper et al. (2008) considered the 
potential of a spherical stellar bulge and an analytical ex- 
pression for disc potential, but no dark matter halo, and an 
ISM that is stratified in z-direction with an effective sound 
speed that is ~ 5 times the normal gas sound speed. Suchkov 
et al. (1994) considered the potential of a spherical bulge and 
a dark matter halo and an initial ISM that is spherically 
stratified. Fragile et al. (2004) considered a spherical halo 
and a z-stratified ISM. However, in a recent simulation of 
outflows driven by supernovae from disc galaxies, Dubois & 
Teyssier (2008) found that the outflowing gas has to contend 
with infalling material from halo, which inhibits the outflow 
for a few Gyr. Fujita et al. (2004) also studied outflows from 
pre-formed disc galaxies in the presence of a cosmological 
infall of matter. 

We choose a z-stratified gas to fill the simulation box, 
with a scale height of 100 pc. For the M2 and M3 case (of 
exponential disc), we also assume a radial profile for the 
initial gas, with a scale length of 5 kpc. For the M3 case, we 
further assume this gas to rotate with decreasing with 
a scale height of 5 kpc. These values are consistent with 
the observations of Dickey & Lockman (1990) and Savage et 
al. (1997) for the warm neutral gas (T ~ 10 4 K) in Milky 
Way. We note that although the scale height for the warm 
neutral gas in our Galaxy is ~ 400 pc at the solar vicinity, 
this is expected to be smaller in the central region because 
of strong gravity due to bulge. The density of the gas just 
above the disc is assumed to be 0.1 particles /cc (0.025 in 
simulation units). 

Furthermore, the adiabatic index of the gas is 5/3 and 



Simulation 



1.5 

> 



0.5 







Total 




Halo 




- Bulge 




Disc 




— Used in this work 






a ' * •»'■"* ' " 




m ** - ** ' 




mZ 

m: y 





,e , , , , 

0.2 0.4 0.6 0.8 1 



r(10kpc) 

Figure 4. Rotation curves corresponding to the gravitational 
fields of an exponential disc, bulge and halo are shown here in 
the units of v Te f[= 200 km s _1 ], along with the total rotation 
curve. The approximation used in our simulation is shown by 
thick red line. 

the gas is assumed initially to be at the same temperature 
corresponding to an initial sound speed c 3 (ini) — 0.1v re f, a 
value which is consistent with the values in our Galaxy for 
the warm ionized gas with sound speed ~ 18 km s~ . 

Our computation domain is rd x rd in the r — z plane, 
with a resolution 512 x 512 cells. The size of individual 
computational cell is ~ 20 pc. We have imposed reflective 
boundary condition around the axis and zero rotational ve- 
locity on the axis. Continuous boundary conditions are im- 
posed at r — rd and z — rd- The lower boundary is slightly 
above the galactic disc with an offset zo = 0.01. We impose 
fixed boundary condition at lower z boundary. The velocity 
of the injected matter is v z (r,zo) = vo = 10 v r e/, and its 
density is given by, 

p{r,z ) = Pzo , forUD (22) 

= P ZQ CXP ' ^ ED ■ 

The density of the injected matter at the base p Z(l = 0.025 
(corresponding to 0.1 protons per cc). 

For the case of exponential disc with rotation (M3), we 
assume for the injected matter to have an angular momen- 
tum corresponding to an equilibrium rotation profile. We 
show in Figure [4] the rotation curves at z = for all compo- 
nents (disc, bulge and halo) separately and the total rotation 
curve. We use the following approximation (shown by thick 
red line in Figure [4| which matches the total rotation curve, 

v (r,z o ) = 1.6 v c [1 - exp(-r/0.15r d )] . (23) 

We assume a bulge of mass Mb = 0.1Af re / and radius 
rt — 0.2L re f. The scale radius for NFW halo (R s ) is de- 
termined for a halo mass Mh = 20Md, as prescribed by 
MMW98. The corresponding disc scale radius is found to be 
r s ~ 5.8 kpc, again using MMW98 prescriptions. There- 
fore we set the disc scale length for the ED case to be 

r s ~ 0.58L r e/. 

The above initial conditions have been chosen to satisfy 
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Table 1. Models. 



Model name 


r 




Disc type 


Mi 


2.0 


0.0 


UD 


M 2 


2.0 


0.0 


ED 


M 3 


2.0 


1.0 


ED 



the following requirements in order to sustain a radiatively 
driven wind as simulated here. 

(i) The strong coupling between dust grains and gas par- 
ticles require that there are of order ~ md/m p number of col- 
lisions between protons and dust grains of mass rrid ~ 10~ 14 
g, for size a ~ 0.1 fim with density ~ 3g cm -3 . To ensure suf- 
ficient number of collisions, the number density of gas parti- 
cles should be n > ^-^-r 1 10~ 3 cm" 3 , for L ref = 10 

kpc. 

(ii) The time scale for radiative cooling of the g 
sumed to be at T ~ 10 4 K, is t cool ~ where A ~ 10~ 23 
erg cm 3 s _1 (Sutherland & Dopita 1993; Table 6) for solar 
metallicity. The typical density filling up the wind cone in 
the realistic case (M 3 ) is ~ 10 _3 -10- 4 cm 3 which gives 
tcooi ~ 8-80 Myr and the dynamical time scale of the wind 
is t re f ~ 50 Myr. Hence radiative cooling is marginally im- 
portant and we will address the issue of radiative cooling in 
a future paper. 

(iii) Radiative transfer effects are negligible since the 
total opacity along a vertical column of length L re f is 
K(nm p )L re f ~ 0.003, for n ~ 10 -3 cm -3 and k ~ 100 cm 2 
g" 1 . 

(iv) The mediation of the radiation force by dust grains 
also implies that the gas cannot be too hot for the dust 
grains to be sputtered. The sputtering radius of grains em- 
bedded in even in a hot gas of temperature T~ 10 5 K is 
~ 0.05(n/0.1 /cc) pm in a time scale of 100 Myr (Tielens et 
al. 1994), and this effect is not important for the tempera- 
ture and density considered here. 



3.2 Simulation set up 

We present 3 models with parameters listed in the Table [T] 
The initial condition for all the models are described in §3.1. 
The boundary condition is essentially same, except that the 
mass flux into the computational domain from the lower z 
boundary depends on the type of disc. As has been men- 
tioned in section 3.1, we keep the velocity of injected matter 
very low, v z (r,zo) = v z (ini) = 10 _5 w re /, so that it does not 
affect the dynamics. The three models have been constructed 
by a combination of different values of three parameters To , 
and the distribution of the density in the disc. Model M3 
has been run for different values of Fo, to ascertain the effect 
of radiation. 

4 RESULTS 

In Figure [5] we present the model Mi for a constant surface 
density disc (UD). The density contour and the velocity vec- 
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Figure 5. Mi : Logarithmic density contours for radiation driven 
wind from UD for four snapshots running up to t = 98 Myr, with 
velocity vectors shown with arrows. Densities are colour-coded 
according to the computational unit of density, 6.7 X 10 -24 g 
cm -3 ~ 4m„ cm -3 . 
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Figure 6. M2 : Logarithmic density contours for radiation driven 
wind from ED for four snapshots running up to t = 98 Myr, with 
velocity vectors shown with arrows. 



tors for the wind are shown in four snapshots in Figure ([5]) 
upto a time t = 98 Myr (corresponding to t — 2 in compu- 
tational time units). There are a few aspects of the gaseous 
flow that we should note here. Firstly, the disc and the out- 
flowing gas in this case has no rotation (v$ — 0). In the 
absence of the centrifugal force due to rotation which might 
have reduced the radial gravitational force, there is a net 
radial force driving the gas inward. At the same time, the 
radiation force, here characterized by Fo = 2, propels the 
gas upward (the radial component of radiation being weak). 
The net result after a few Myr is that the gas in the region 
near the pole moves in the positive z direction, and there is 
a density enhancement inside a cone around the pole, away 
from which the density and velocities decrease. 

Also, because of the strong gravity of the bulge, the 
gas tends to get trapped inside the bulge region, and even 
the gas at larger r tends to get dragged towards the axis. 
This region puffs due to accumulation of matter. Ultimately 
the radiative force drives matter outwards in the form of a 
plume. 

Next, we change the disc mass distribution and simulate 
the case of wind driven out of an exponential disc (ED). 
We show the results in Figure [6] Since both gravity and 
radiation forces in this case of exponential disc are quasi- 
spherical in nature, therefore in the final snapshot the flow 
appears to follow almost radial streamlines. Although in the 
vicinity of the disc, the injected matter still falls towards the 
axis, but this is not seen at large height as was seen in the 
previous case of Mi . This makes the wind cone of rising gas 
more diverging than in the case of UD (Mi). 



4.1 Rotating wind from exponential disc 

The direction of the fluid flow in Mi and M2 is by and large 
towards the axis, and this flow is mitigated in the presence of 
rotation in the disc and injected gas. In the next model M3, 
we consider rotating matter being injected into the compu- 
tational domain and which follows a distribution given 
by Eq. (|23p . This is reasonable to assume since the disc from 
which the wind is supposed to blow, is itself rotating. In M3, 
we simulate rotating gas being injected above a ED and be- 
ing driven by a radiation force of To = 2. We present nine 
snapshots of the M3 case in Figure [7] 

The first six snapshots of Figure[7]show the essential dy- 
namics of the outflowing gas. The fast rotating matter from 
the outer disc is driven outward because the radial gravity 
component is balanced by rotation. Near the central region, 
rotation is small and also the radial force components are 
small. Therefore the gas is mostly driven vertically. The in- 
jected gas reaches a vertical height of ~ 5 kpc in a time scale 
of ~ 37 Myr (t=0.75). The flow reaches a steady state after 
~ 60 Myr (t=1.25). In the steady state we find a rotating 
and mildly divergent wind. 

We show the azimuthal velocity contours in Figure [8] in 
colour for the fully developed wind (last snapshot in M3), 
and superpose on it the contour lines of p. The density con- 
tours clearly show a conical structure for outflowing gas. 
The rotation speed of the gas peaks at the periphery of the 
cone, and is of order ~ 50-100 km s . Compared to the 
disc rotation speed, the rotation speed of the wind region is 
somewhat smaller. In other words, we find the wind mostly 
consisting of low-angular momentum gas lifted from the disc. 

We plot the velocity of gas close to the axis in Figure 
[9] for different times in this model (M3), using v(0, z) ~ 
v z (0 + ,z). The velocity profile in the snapshots at earlier 
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Figure 7. M3: Contours of log 10 (p) and v-field of radiation driven wind with To = 2.0 from an ED. t = 2 corresponds to 98 Myr. 



v^(200 km/s) 
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r(10kpc) z (10 kpc) 

Figure 8. The rotation velocity for the case M3 at a time of Figure 9. The axial velocity v z (0 + ,z) with 2 at different time 

98 Myr is shown in colours. Contour lines of logio(p) are plotted steps for the model M3. t = 2.0 corresponds to a time of 98Myr. 
over it. 
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Figure 10. The axial velocity n z (0 + , lOfcpc) in simulation units 
v re f = 200 km s — 1 with Tq, at a time t ~ 10 2 Myr. 



time fluctuates at different height, but becomes steady after 
t 1.5, as does the density profile. 

We have run this particular case of ED with rotation 
(model M3) for different values of To. In order to illustrate 
the results of these runs, we plot the z-component of velocity 
(fz(0 + ,10 kpc)) at 10 kpc and at simulation time, t — 2 
as a function of To in Figure 1101 We find that significant 
wind velocities are obtained for To > 1.5 and wind velocities 
appear to rise linearly with To after this critical value is 
acheived. Sharma et al. (2011) found this critical value to 
be To ~ 2 for a constant density disc and wind launched 
above the bulge. For the realistic case of an exponential disc, 
we find in the present simulation the critical value to be 
somewhat smaller than but close to the analytical result. 
The important point is that the critical Fo is not unity. 
This is because the parameter To is not a true Eddington 
parameter since it is defined in terms of disc gravity and 
radiation, whereas halo and bulge also contribute to gravity. 



5 DISCUSSIONS 

Our simulation differs from earlier works (e.g. Suchkov et.al. 
1994) mainly in that we specifically target warm outflows 
and the driving force is radiation pressure. Most of the pre- 
vious simulations of galactic wind have used energy injected 
from supernovae blasts as a driving force. However, with the 
ideas presented in Murray et al. (2005), which worked out 
the case of radiation pressure in a spherical symmetric set- 
up, it beomes important to study the physics of this model 
in an axisymmetric set up, as has been done analytically by 
Sharma et al. (2011) (see also, Zhang & Thompson 2010). 
Also we have tried to capture all features of a typical disc 
galaxy like a bulge and a dark matter halo, and a rotating 
disc. Recent analytical works (Sharma & Nath 2011) and 
simulations (Hopkins et al. 2011) have shown that outflows 



from massive galaxies {Mhalo ^ 10 12 Mq) have different 
characteristics than those from low mass galaxies. Outflows 
from massive galaxies are mostly driven by radiation pres- 
sure and the fraction of cold gas in the halos of massive 
galaxies is large (van de Voort & Schaye 2011). Our simula- 
tions presented here addresses these outflows in particular. 

We have parameterized our simulation runs with the 
disc Eddington factor To, and it is important to know the 
corresponding luminosity for a typical disc galaxy, or the 
equivalent star formation rate. For a typical opacity of a 
dust and gas mixture [k ~ 200 cm 2 g _1 ) (Draine 2011), the 
correspondig mass-to-light ratio requirement for Fo > 1.5 is 
that M/L < 0.03. Sharma et al. (2011) showed that for the 
case of an instantaneous star formation, To > 2 is possible 
for an initial period of ~ 10 Myr after the starburst. However 
for a continuous star formation, which is more realistic for 
disc galaxies, Sharma & Nath (2011) found that only ultra 
luminous infrared galaxies (ULIGs), with star formation rate 
larger than ~ 100 Mq yr _1 and which are also massive, 
are suitable candidates for such large values of Fo, and for 
radiatively driven winds. 

The results presented in the previous sections show that 
the outflowing gas within the central region of a few kpc 
tends to stay close to the pole, and does not move out- 
wards because of its low angular momentum. This makes 
the outflow somewhat collimated. Although outflows driven 
by SN heated hot wind also produces a conical structure 
(e.g., Fragile et al. 2004) emanating from a breakout point 
of the SN remnants, there is a qualitative difference between 
this case and that of radiatively driven winds as presented 
in our simulations. While it is the pressure of the hot gas 
that expands gradually as it comes out of a stratified at- 
mosphere, in the case of a radiation driven wind, it is the 
combination of mostly the lack of rotation and almost verti- 
cal radiation driving force in the central region that produce 
the collimation effect. 

We also note that the conical structure of rotation in 
the outflowing gas is similar to the case of outflow in M82 
(Greve 2004), where one observes a diverging and rotating 
periphery of conical outflow. 

We have not considered radiative cooling in our simu- 
lations, since for typical density in the wind the radiative 
cooling time is shorter or comparable than the dynamical 
time. However, there are regions of higher density close to 
the base and radiative cooling can be important there. We 
will address this point in a future paper. 

From our results of the exponential and rotating disc 
model, we find the wind comprising of low-angular mo- 
mentum gas lifted from the disc. It is interesting to note 
that recent simulations of supernovae driven winds have also 
claimed a similar result (Governato et al. 2010). Such loss 
of low angular momentum gas from the disc may have im- 
portant implication for the formation and evolution of the 
bulge, since the bulge population is deficient in stars with 
low specific angular momentum. Binney, Gerhardt & Silk 
(2001) have speculated that outflows from disc that pref- 
erentially removes low angular momentum material may re- 
solve some discrepancies between observed properties of disc 
and results of numerical simulations. 

As a caveat, we should finally note that the scope and 
predictions of our simulation is limited by the simple model 
of disc radiation adoped here. In reality, radiation from disks 
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is likely to be confined in the vicinity of star clusters, and not 
spread throughout the disk as we have assumed here. This 
is likely to increase the efficacy of radiation pressure, but 
which is not possible within the scope of an axisymmetric 
simulation. 



6 SUMMARY 

We have presented the results of hydrodynamical (Eulerian) 
simulations of radiation driven winds from disc galaxies. 
After studying the cases of winds from a constant surface 
density disc and exponential disc without rotation, we have 
studied a rotating outflow originating from an exponential 
disc with rotation. We find that the outflow speed increases 
rapidly with the disc Eddington parameter To — k//(2cGE) 
for To ^ 1.5, consistent with theoretical expectations. The 
density structure of the outflow has a conical appearance, 
and most of the ouflowing gas consists of low angular mo- 
mentum gas. 

We thank Yuri Shchekinov for constructive comments 
and critical reading of the manuscript. IC acknowledges the 
hospitality of the Astronomy and Astrophysics Group of Ra- 
man Research Institute, where the present work was con- 
ceived. DR was supported by National Research Foundation 
of Korea through grant 2007-0093860. 



REFERENCES 

Bianchi, S., Ferrara, A. 2005, MNRAS, 358, 379 
Binney, J., Gerhard, O., Silk, J. 2001 MNRAS, 321, 471 
Burke, A. J. 1968, MNRAS, 140, 241 
Chattopadhyay, I. 2005, MNRAS, 356, 145 
Chevalier, R. A., Clegg, A. W., 1985, Nature, 317, 44 
Chiao, R. Y., Wickramasinghe, N. C, 1972, MNRAS, 159, 
361 

Cooper, J. L., Bicknell, G. V., Sutherland, R. S., Bland- 
Hawthorn, J. 2008, ApJ, 674, 157 

Dave, R., Finlator, K., Oppenheimer, B. D. 2006, MNRAS, 
370, 273 

Davies, J. I., Alton, P., Bianchi, S., Trewhella, M. 1998, 

MNRAS, 300, 1006 
Dekel, A., Silk, J. 1986 ApJ, 303, 39 
Dickey, J. M., Lockman, F. J. 1990, ARA&A, 28, 215 
Draine, B. T. 2011, ApJ, 732, 100 
Dubois, Y., Teyssier, R. 2008 A&A, 477, 79 
Fragile, P. C, Murray, S. D., Lin, D. N. C. 2004, ApJ, 617, 

1077 

Fujita, A., Mac Low, M., Ferrara, A., Meiksin, A. 2004, 
ApJ, 613, 159 

Fujita, A., Martin, C. L., Mac Low, M., Kimberly, C. B., 

Weaver, R. 2009, ApJ, 698, 693 
Governato et. al. 2010, Nature, 463, 203 
Greve, A. 2004, A&A, 416, 67 
Harten, A. 1983, J. Comput. Phys., 49, 357 
Heckman, T. M. 2002, Extragalactic Gas at Low Redshift 

(ASP Conf. Ser. 254), ed. J. S. Mulchaey & J. T. Stocke 

(San Fransisco, CA: ASP), 292 
Hopkins, P. F., Quataert, E., Murray, N. 2011, preprint 

(|arXiv:1110.4638|l 
Johnson, H. E., Axford, J. I. 1971 ApJ, 165, 381 



Kang, H., Ostriker, J. P., Cen R., Ryu, D., Hernquist, L., 
Evrard, A. E., Bryan G. L., & Norman, M. L. 1994, ApJ, 
430, 83 

Kohji, T., Ikeuchi, S. 1988, ApJ, 330, 695 
Larson, R. B. 1974, MNRAS, 169, 229 
Lynds C, Sandage A. 1963 ApJ, 137, 1005 
Mac Low, M., Ferrara, A. 1999, ApJ, 513, 142 
Martin, C. L. 2005, ApJ, 621, 227 
Mathews, W. G., Baker, R. G. 1971, ApJ, 170, 241 
Mo, H. J., Mao, S., White, S. D. M. 1998, MNRAS, 295, 
319 (MMW98) 

Molteni, D., Ryu, D., Chakrabarti, S. K. 1996, ApJ, 470, 
460 

Murray, N., Quataert, Q. & Thompson, T. A. 2005, ApJ, 
618, 569 

Murray, N., Menard, B., Thompson, T. A. 2011, ApJ, 735, 
66 

Nath, B. B. & Silk, J. 2009, MNRAS, 396, L90 
Navarro, J. F., Frenk, C. S., White, S. D. M. 1997, ApJ, 
490, 493 

Oppenheimer, B. D., Dave, R. 2006, MNRAS, 373, 1265 
Rupke, D. S., Veilleux, S., Sanders, D. B. 2005, ApJS, 160, 
115 

Rupke, D. S., Veilleux, S., Sanders, D. B. 2002, ApJ, 570, 
588 

Ryu, D., Ostriker, J. P., Kang, H., Cen, R., 1993, ApJ, 414, 
1 

Savage, B. D., Sembach, K. R., Lu, L. 1997, AJ, 113, 2158 

Seaquist, E. R., Clark, J. 2001, ApJ, 552, 133 

Shapley, A. E., Steidel, C. C, Pettini, M. & Adelberger K. 

L. 2003, ApJ, 588, 65 
Sharma, M., Nath, B. B., Shchekinov, Y. 2011, ApJ, 736, 

L27 

Sharma, M., Nath, B. B. 2011, preprint (|arXiv:1112.3447|> 
Sofue, Y., Reuter, H. P., Krause, M., Wielebinski, R., 

Nakai, N. 1992, ApJ, 395, 126 
Strickland, D. K., Stevens, I. R. 2000, ApJ, 314, 511 
Strickland, D.K., Heckman, T.M., Colbert, E.J.M., 

Hoopes, C.G., Weaver, K.A. 2004 ApJS, 151, 153 
Suchkov, A. A., Berman, V. G., Heckman, T. M., Balsara, 

D. S. 1994, ApJ, 430, 511 
Suchkov, A. A., Berman, V. G., Heckman, T. M., Balsara, 

D. S. 1996, ApJ, 463, 528 
Sutherland, R. M., Dopita, M. A. 1993, ApJS, 88, 253 
Tielens, A. G. G. M., McKee, C. F., Seab, C. C, Hollen- 

bach, D. J. 1994, ApJ, 431, 321 
Tomisaka, K., Bregman, J. N. 1993, PAS J, 45, 513 
van de Voort, F., Schaye, J. 2011, preprint 

(|arXiv:1111.5039|> 
Veilleux S., Cecil G., Bland-Hawthorn J., 2005, ARA&A, 

43, 769 

Veilleux S., Rupke, D. S. N., Swaters, R. 2009, ApJ, 700, 
L149 

Walter, F., Weiss, A., Scoville, N. 2002, ApJ, 580, L21 
Wesmoquette et al., 2009, ApJ, 696, 192 
Zhang, D., Thompson, T. A. 2010, 
preprint (|arXiv: 1005.4691 P 



10 /. Chattopadhyay, M. Sharma, B. B. Nath and D. Ryu 



7 APPENDIX 

Consider a razor thin disc in rcj> plane as illustrated in the 
fig. Now our task is to calculate the force components at any 
arbitrary point above the disc. Let us consider an annulus 
of the disc between rl and rl + drl. Area of the element at 
point P(r7, (f>/,0) is rldrld(f>l. Also take a field point Q(r,0,z) 
above the disc plane. Azimuthal coordinate of Q is taken 
to be zero for simplicity as we know that azimuthal force 
components are zero due to symmetry. Let QN and QM be 
perpendiculars from Q on the x and the z axis, respectively. 
So we can write, 



PN 
PQ 2 



(r — r/cos(f)/) + (rlsirufil) 



PN 2 + z 2 



2 , 2 , ,2 

r + z + rl 



2rrlcos(j)l 



PN 

sinZPQN = — 



(24) 



The gravitational force due to the small area element at P 
is given by 

G dm PQ , 



dF 



(PQ) 3 



-n ; 



dm = ridrid<piT,(ri) 



(25) 



Here E(r/) is the surface density of the disc. Now the z com- 
ponent of this force is 



zGY.(rl) r!drtd(j)t 



[r 2 + z 2 + rl 2 — 2rr/cos(p/] a i 12 



(26) 



dF g , z = \dF s \ — 

To calculate the radial component, let PS be the perpendic- 
ular from P on the x-axis. Then, we have sinZSPN = SN/PN 
= (r-r/cos0/)/PN. Component of the force along the direc- 
tion of PN is 



dF g , PN = \dF g \sinZPQN = 
So the radial component is 

SN PN 

dF g , r = d F 3 .p N sinASPN = \dF s \ — — 



(27) 



(r - rtcostf)/) GE(r/) rldrldfyl 
2rrlcos(f>i\ i / 2 



(28) 




Figure 11. Schematic diagram for the calculation of gravita- 
tional force due to disc in the icy-plane. We consider an annulus 
in the disc and an element of area around the point P (rl, <j>, 0) in 
this annulus is considered here in order to compute the force at a 
point Q (r, 0, z) whose azimuthal coordinate = 0. The point S 
(rl cos 0,0,0) is the foot of the perpendicular drawn from P on the 
x-axis. The point Si (rl, 0, z) is at the intersection of the vertical 
from S (along z-axis) and the line parallel to x-axis at height z. 
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